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Abstract 

We  apply  an  inverse  problem  formulation  to  determine  characteristics  of  a  defect  from  a  perturbed 
electromagnetic  interrogating  signal.  A  defect  (gap)  inside  of  a  dielectric  material  causes  a  disruption, 
via  reflections  and  refractions  at  the  material  interfaces,  of  the  windowed  interrogating  signal.  We  model 
the  electromagnetic  waves  inside  the  material  with  Maxwell’s  equations.  This  leads  to  a  non-standard, 
nonlinear  optimization  problem  for  the  dimensions  and  location  of  the  defect.  Using  simulations  as 
forward  solves,  we  employ  a  Newton-based,  iterative  optimization  scheme  to  a  novel  modified  least- 
squares  objective  function.  Numerical  results  are  given  in  tables  and  plots,  standard  errors  are  calculated, 
and  computational  issues  are  addressed. 


1  Introduction 

The  problem  we  consider  is  that  of  detecting  a  gap  inside  of  a  dielectric  material  using  high  frequency 
electromagnetic  interrogation.  The  idea  is  to  observe  the  reflected  and/or  transmitted  signals  and  use 
the  data  to  solve  an  inverse  problem  to  determine  certain  characteristics  of  the  gap,  e.g.,  location  and/or 
width.  Possible  applications  of  this  procedure  include  quality  assurance  in  fabrication  of  critical  dielectric 
materials,  or  damage  detection  in  aging  materials  for  safety  concerns.  Further  applications  of  electromagnetic 
interrogation,  as  well  as  additional  problem  formulations  and  solutions,  can  be  found  in  [2]. 

The  particular  motivation  for  this  research  is  the  detection  of  defects  in  the  insulating  foam  on  the  fuel 
tanks  of  the  space  shuttles  in  order  to  help  eliminate  the  separation  (delamination)  of  foam  during  shuttle 
ascent.  To  this  end  we  address  the  problem  of  detecting  a  single  gap  formed  between  a  dielectric  medium 
and  a  supra-conducting  backing  representing  the  foam  and  the  metallic  tank,  respectively.  However,  first  we 
develop  our  methodology  on  the  slightly  simpler  problem  of  a  gap  formed  in  the  interior  of  the  foam  (void) , 
where  for  simplicity,  we  ignore  the  reflections  from  the  back  boundary  (i.e.,  we  impose  absorbing  boundary 
conditions  instead) .  We  also  allow  for  the  possibility  in  this  case  that  the  foam  has  been  removed  for  testing, 
and  therefore  we  are  able  to  place  sensors  both  on  the  front  and  back  sides  of  the  foam. 

To  be  applicable  to  real  world  problems  we  must  eventually  be  able  to  solve  these  inverse  problems  with 
length  scales  on  the  order  of  20cm  for  the  thickness  of  the  foam,  .2 mm  for  the  width  of  the  gap,  and  a 
wavelength  of  about  3 mm.  This  wavelength  corresponds  to  a  frequency  of  100 GHz,  which  is  the  lower  end 
of  the  terahertz  frequency  range  (.1  ~  10  x  10 12Hz).  The  rationale  for  using  this  choice  of  frequency  is 
that  higher  frequencies  are  significantly  attenuated  in  the  materials  which  we  are  interested  in  interrogating. 
Lower  frequencies  (larger  wavelengths)  have  less  resolution  in  detecting  small  gaps,  and  are  less  capable  of 
sharply  distinguishing  between  air  and  foam  which  may  have  a  high  air  content. 

First  we  simplify  the  problem  by  considering  a  linearly  polarized,  pulsed  interrogating  signal  which 
reduces  the  problem  to  one  spatial  dimension.  We  then  define  an  inverse  problem  for  determining  the  gap’s 
dimensions.  We  assume  that  we  have  data  from  sensors,  located  in  front  of  and/or  behind  the  material,  that 
record  the  electromagnetic  signal  after  it  is  reflected  from  (or  passes  through)  the  material  interfaces.  We 
compute  simulated  signals  with  approximations  to  the  gap’s  characteristics  and  apply  an  optimization  routine 
to  a  modified  least  squares  error  between  this  simulated  signal  and  the  given  data.  In  our  computations  we 
use  Maxwell’s  equations  and  a  Debye  polarization  equation  to  model  the  signal,  and  solve  these  equations 
using  a  Finite  Element  method  in  space  and  Finite  Difference  methods  in  time.  Thus  the  optimization 
routine  finds  those  gap  characteristics  which  generate  a  simulated  signal  that  most  closely  matches  the  given 
data.  In  this  sense  we  determined  an  estimate  to  the  “true”  gap  characteristics. 
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Figure  1: 


Problem  1 :  Dielectric  slab  with  a  material  gap  in  the  interior.  Possible  sensors  in  front  and  behind. 


In  Section  2  we  define  the  equations  that  we  have  chosen  in  order  to  model  the  electromagnetic  waves 
inside  the  material.  We  further  distinguish  between  two  problem  types  that  we  will  address  (denoted  as 
Problem  1  and  Problem  2).  Section  3  contains  the  details  of  our  numerical  methods  for  the  simulations.  We 
introduce  the  inverse  problem  formulation  for  Problem  1  in  Section  4,  and  later  improve  upon  it  in  Section 
4.2.  Numerical  results  of  the  inverse  problem  are  displayed  in  Section  4.3. 

In  Section  5  we  begin  addressing  Problem  2.  Similarities  and  differences  between  the  computational 
issues  between  the  two  problems  are  pointed  out.  A  more  sophisticated  optimization  method  is  described  in 
Section  5.3  and  associated  numerical  results  are  given  in  Section  5.5.  In  Sections  5.5.1  and  5.5.2  we  explore 
the  effects  of  adding  random  noise  to  the  data,  both  relative  and  constant  variance.  In  the  latter,  we  compute 
standard  error  estimates. 


2  Problem  Description 


We  interrogate  an  (infinitely  long)  slab  of  homogeneous  nonmagnetic  material  by  a  polarized,  windowed 
signal  (see  [2]  for  details)  in  the  THz  frequency  range  (see  Figure  1).  We  assume  a  wave  normally  incident 
on  a  slab  which  is  located  in  O  =  [21,24]  with  faces  parallel  to  the  x-y  plane  (see  Figure  2).  Note  that  we 
employ  the  “method  of  mappings”  (see  [2])  in  our  computations,  therefore  we  may  assume  0  <  21  <  24  <  1. 
We  denote  the  vacuum  outside  of  the  material  by  fl0-  The  electric  field  is  polarized  to  have  oscillations  in 
the  x-z  plane  only.  Restricting  the  problem  to  one  dimension,  we  can  write  the  electric  and  magnetic  fields, 
E  and  H  respectively,  as  follows 

E(t,  x)  =  iE(t ,  2) 

H(t,x)  =  jH(t,  2), 


so  that  we  are  only  concerned  with  the  scalar  values  E(t,  2)  and  H(t,  2). 
Maxwell’s  equations  then  become: 


dE 

dz 

dH 

dz 


dH 
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Figure  2:  The  domain  of  the  material  slab:  O  =  [z±,  Z4]. 


where  D(t,  z)  is  the  electric  flux  density,  /io  is  the  magnetic  permeability  of  free  space, 
and  Js  is  a  source  current  density  (determined  by  our  interrogating  signal).  We  take 
of  Equation  (la)  with  respect  to  z,  and  the  partial  of  Equation  (lb)  with  respect  to 
terms  in  each,  and  thus  eliminating  the  magnetic  field  H ,  we  have: 

E"  =  fx0  (D  +  aE  +  js ^  , 

(where  '  denotes  z  derivatives  and '  denotes  time  derivatives) . 

Note  that  we  have  neglected  magnetic  effects  and  we  have  let  the  total  current  density  be  J  =  Jc  +  Js, 
where  Jc  =  aE  is  the  conduction  current  density  given  by  Ohm’s  law  in  the  material. 

For  our  source  current,  Js,  we  want  to  simulate  a  windowed  pulse,  i.e.,  a  pulse  that  is  allowed  to  oscillate 
for  one  full  period  and  then  is  truncated.  Further,  we  want  the  pulse  to  originate  only  at  z  =  0,  simulating 
an  infinite  antenna  at  this  location.  Thus  we  define 

Js{t,z)  =  S(z)sin(wt)I[o!tf](t), 

where  ui  is  the  frequency  of  the  pulse,  tf  =  2'k/lo  is  fixed,  I[o.tf]  (t)  represents  an  indicator  function  which  is 
1  when  0  <  t  <  tf  and  zero  otherwise,  and  d(z)  is  the  Dirac  delta  distribution. 

Remark  1  Computationally,  having  a  windowed  signal  introduces  discontinuities  in  the  first  derivatives 
which  are  not  only  problematic  in  the  numerical  simulations  (producing  spurious  oscillations),  but  are  also 
essentially  non-physical.  Therefore  in  our  implementation  we  actually  multiply  the  sine  function  by  an  expo¬ 
nential  function  (see  [3]  for  details)  rather  than  the  traditional  indicator  function.  However,  for  notational 
consistency  we  will  continue  to  denote  this  function  as  /[0  t/](t). 

The  electric  flux  density  inside  the  material,  given  by  D  =  eoCoo E  +  P,  is  dependent  on  the  polarization, 
P.  Note  that  eo  is  the  permittivity  of  free  space  and  e^  is  the  relative  permittivity  in  the  limit  of  high 
frequencies.  For  computational  testing  we  assume  for  this  presentation  that  the  media  is  Debye  and  thus  we 
use  the  following  polarization  model  inside  O: 

TP  +  P  =  e0(es  -  Coo )E, 

where  es  is  a  static  relative  permittivity  and  r  is  a  relaxation  time.  We  also  assume  P(0,z)  =  0.  Note  that 
in  the  vacuum  outside  of  f2,  P  =  0. 


a  is  the  conductivity, 
the  partial  derivative 

r\  2  T  T 

t.  Equating  the 
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In  order  to  represent  D  in  the  entire  domain,  we  use  the  indicator  function  Iq  which  is  1  inside  17  and 
zero  otherwise.  Thus 

D  =  eo  E  +  eo(£oo  —  1  )Iq.E  +  IqP. 

In  order  to  have  a  finite  computational  domain,  we  impose  absorbing  boundary  conditions  at  2  =  0  and 
2  =  1,  which  are  modeled  as 


E-cE' 
E  +  cE' 


=  0 


J  z=0 


=  0. 


With  these  boundary  conditions,  any  boundary  incident  signal  passes  out  of  the  computational  domain,  and 
does  not  return,  i.e.,  we  force  it  to  be  absorbed  by  the  boundary.  Also  we  assume  zero  initial  conditions,  i.e., 


E(0,z)  =  0 
E(0,z)  =  0. 


(2) 


Thus  our  entire  system  can  be  written 

MoCo(l  +  (eoo  ^  1  )In)E  +  po IqP  +  novInE  -  E"  =  ~n0js  in  17  U(l0 

tP  +  P  =  e0(es  -  eoo)E  in  O 
[E  -  cE']z-o  =  0 
[E  +  cE’}z=1  =  0 


(3) 


with  (2)  and 

Js(t,z)  =  6(z)sin(ujt)I[0ttf](t).  (4) 

Classical  solutions  to  (3)  should  not  be  expected  due  to  the  windowed  interrogating  signal  and  the 
discontinuous  dielectric  parameters  across  interfaces.  For  this  reason,  and  also  to  enable  the  application 
of  the  Finite  Element  Method,  we  prefer  to  convert  (3)  to  weak  form  using  spaces  H  =  £2(0,1)  and 
V  =  H1( 0, 1).  Substituting  er  =  (1  +  (e^  —  1)/q)  and  =  es  —  results  in  the  following  weak  system 

(n0e0erE,  (j>)  +  {no IqP,  </>)  +  {hquIqE,  (/>)  +  (E1 ,  <//)  -  -E(t,  1)0(1)  +  -E(t,  0)0(0)  =  -(/x0  js,  0) 

c  c  (5) 

(tP,  0)n  +  (P,  0)o  =  (e0 edE,  0)n 


with  (2)  and  (4).  Note  that  (•,  •)  is  modified  from  the  traditional  L2  inner  product  due  to  the  aforementioned 
use  of  the  “method  of  mappings”  (see  [2]).  Existence  and  uniqueness  of  systems  of  this  type  are  treated  in 
[21- 

In  this  formulation  we  have  initially  assumed  a  single  slab  of  a  dielectric  contained  in  17  =  [21,24]. 
Thus  Iq  =  1  if  21  <  2  <  24,  and  zero  otherwise.  We  now  introduce  a  gap  consisting  of  a  vacuum  in 
the  interior  of  the  material  as  depicted  in  Figure  3.  If  the  gap  is  located  in  (22,23)  then  we  redefine 
O  =  {2I21  <  2  <  22  or  23  <  2  <  24}.  We  will  refer  to  this  formulation  as  Problem  1  (recall  Figure  1). 
Note  that  it  is  not  necessary  to  enforce  additional  conditions  at  the  interfaces  as  they  are  natural  interface 
conditions  and  are  implied  within  the  weak  form  of  the  system.  Later  we  will  discuss  a  second  problem 
formulation,  Problem  2,  where  the  gap  is  between  the  dielectric  slab  and  a  metallic  (supra-conducting) 
backing,  as  shown  in  Figure  4.  This  will  require  slightly  different  boundary  conditions  (reflecting  instead 
of  absorbing  at  2  =  1,  where  the  metal  backing  begins),  but  otherwise  the  numerical  solution  methods  and 
analysis  are  the  same. 


3  Numerical  Solution 

3.1  Finite  Elements 

We  apply  a  Finite  Element  method  using  standard  linear  one  dimensional  basis  elements  to  discretize  the 
model  in  space.  Let  N  be  the  number  of  intervals  in  the  discretization  of  2,  and  h  =  1  /IV,  then  the  Finite 
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Figure  3:  The  domain  of  the  material  slab  with  an  interior  gap  between  zi  and  23:  O  =  {z\z±  <  z  < 

Z2  or  23  <  z  <  24}  . 
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Figure  4:  Problem  2:  Dielectric  slab  and  metallic  backing  with  a  gap  between.  Possible  sensors  only  in  front. 


Element  discretization  has  an  order  of  accuracy  of  0(h2).  For  implementation  we  scale  time  by  t  =  ct  and 
the  polarization  by  P  =  P/e 0  for  convenience.  The  resulting  system  of  ordinary  differential  equations  after 
the  spatial  discretization  is  the  semi-discrete  form 


erMe  +  Mnp  +  ( r/oaMn  +  D  +  B)e  +  Ke  =  rj^J  (6a) 

Mnp  +  \Mnp  =  ed\ Mne,  (6b) 
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where  A  =  A ,  and  770  =  \J po  / eo .  Also  e  and  p  are  vectors  representing  the  approximate  values  of  E  and  P, 
respectively,  at  the  nodes  27  =  ih.  The  mass  matrix  M  has  entries 


Mij  =  ~  /  4>i(t>jdz, 

Jo 


where  are  the  basis  functions  (Mn  is  the  mass  matrix  integrated  only  over  Q),  while  the  stability 

matrix  K  has  entries 


Kij  =  ■=  [  fiifijdz. 

Jo 


The  matrices  D  and  B  result  from  the  boundary  conditions  where 


Cl,!  =1 
Bn+i,n+i  =  1 


and  all  other  entries  are  zero.  Finally,  J  is  defined  as 


Ji  —  • —  /  Jsfiidz. 

Jo 

Note  that  by  differentiating  (6b)  we  can  substitute  into  (6a)  and  obtain  an  equation  only  dependent 
explicitly  on  P  (two  substitutions  are  required  to  eliminate  P  and  P ): 

erMe  +  (% aMn  +  D  +  B  +  edXMn)e  +  {K  -  edX 2Mn)e  +  X 2Mnp  =  p0J. 

Using  shorthand  we  can  write  our  entire  coupled  system  as 

M{e  +  A^e  +  A^e  +  X2p  =  7/0  J  (7a) 

p  +  Xp  =  edXMne,  (7b) 

where  p  =  M‘lp.  It  is  important  to  mention  that  each  matrix  is  tridiagonal  due  to  the  choice  of  the  linear 
finite  elements. 


3.2  Finite  Differences 

In  order  to  solve  the  semi-discrete  form  of  our  equations  we  consider  two  distinct  finite  difference  methods. 
In  the  first  method  we  convert  the  coupled  second  order  system  of  equations  into  one  larger  first  order  system 
and  simply  apply  a  theta  method  (unless  otherwise  stated,  we  use  9  =  t).  In  the  second  method  we  solve 
first  for  the  polarization  with  a  forward  differencing  scheme  using  the  initial  conditions  and  then  use  that 
to  update  a  second  order  central  difference  scheme  for  the  magnitude  of  the  electromagnetic  field.  We  then 
continue  this  process  iteratively,  alternating  between  solving  for  P  and  for  E. 

Both  methods  are  second  order  in  time  and  space  for  appropriately  smooth  data  (and  with  At  =  O(h)). 
We  have  compared  the  errors  and  the  run  times  of  the  two  methods  for  several  smooth  test  problems  and 
have  determined  that  the  second  method  is  as  accurate,  but  twice  as  fast,  as  the  first  in  all  cases  primarily 
due  to  the  fact  that  the  linear  system  is  of  smaller  dimension  in  the  second  method.  Also,  the  first  method 
incidentally  solves  for  e  in  addition  to  e  and  p,  which  is  superfluous. 

In  our  second  method  we  use  a  second  order  central  difference  scheme  to  solve  (7a).  Thus  we  must  first 
find  an  approximation  to  E(ti,z)  where  t,;  =  iAt.  Note  that  approximating  E  with  its  Taylor  expansion 
around  to  =  0  and  applying  the  initial  conditions  and  ODE,  one  obtains 

At2 

E(ti,z)  «  -  —  p0Js(0,2). 

Our  approach  is  to  first  solve  for  p  using  a  0-method,  and  then  use  that  approximation  to  solve  for  e  at 
the  next  time  step.  Thus,  our  finite  difference  approximation  for  (7b)  is 

Pn+l  =  Pn  +  1  XAtO  e’n+fl  _  Pn)  (®) 
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where  [en]j  =  E(tn,Zj),  [pn\j  =  Mn P(tn,Zj),  Zj  =  jh,  and  en+g  =  Oen  +  (1  —  0)en+i  is  a  weighted  average 
of  en  and  en+\  for  relaxation  to  improve  the  stability  of  the  method.  Once  we  have  pn+\  we  can  solve  for 
e„+2-  Applying  a  second  order  central  difference  method  with  averaging  to  (7a)  gives 

Aien+2  =  ^Sn+i  +  A3  en  +  At2r/o  Jn+\  —  X2At2pn+i.  (9) 

Note  that  in  this  case  A\  is  tridiagonal  and  the  matrix  is  the  same  for  each  time  step,  so  we  may  store  the 
Crout  LU  factorization  and  use  back  substitution  to  solve  the  system  at  each  time  step.  For  tridiagonal 
matrices  the  factorization  and  the  back  substitution  are  both  order  O(N).  See  [3]  for  computational  issues 
encountered  in  implementation,  including  the  use  of  the  method  of  maps  to  allow  for  gap  sizes  (5)  smaller 
than  the  mesh  size  ( h ). 

3.3  Numerical  Simulations 

The  following  figure  depicts  the  numerical  solution  of  the  amplitude  of  the  electric  field  at  various  times 
(Figure  5).  We  considered  a  Debye  medium  with  the  following  parameters: 

78.2, 

5.5, 

1  x  1CT5, 

3.16  x  1(T8, 

2  GHz, 

and  a  material  gap  located  at  [23,  24]  =  [.6,  .61]. 

We  have  used  2  GHz  simply  so  that  our  computational  domain  of  2  €  [0, 1]  would  not  have  to  be  scaled 
for  this  demonstration.  Also,  in  practice,  one  would  not  compute  a  domain  so  much  larger  than  the  material, 
just  as  in  an  experiment  the  sensors  should  be  as  close  to  the  material  as  possible  to  reduce  noise.  We  did  so 
here  merely  so  that  the  full  wavelength  of  the  signal  would  be  visible.  See  [3]  for  figures  showing  the  signal 
recorded  at  receivers  located  at  2  =  0  and  2  =  1. 


4  Problem  1 

We  now  apply  an  optimization  routine  to  the  least  squares  error  between  a  simulated  signal  and  the  given 
data  to  try  to  determine  the  gap  characteristics.  In  particular  we  will  be  trying  to  find  the  depth,  d  :=  22  —  24, 
and  the  width,  <5  :  =  23  —  22,  which  will  produce  a  simulated  signal  most  closely  similar  (in  the  least  squares 
sense)  to  the  data.  Existence  and  continuous  dependence  of  inverse  problems  of  this  type  are  addressed  in 
[2]- 

4.1  Inverse  Problem 

All  of  the  following  are  solved  with  respect  to  a  reference  problem  (Rl)  with  these  parameter  values  (see 
also  Figure  3): 

20  =  0, 21  =  .2, 22  =  .3, 23  =  .5, 24  =  .8, 25  =  1.0 
/  =  4  GHz,  tf  is  one  period 

r  =  8.1  x  10— 12 ,  a  =  lx  10-5,  es  =  80.1,  £00  =  5.5  in  the  material 
N  =  1024,  Nt  =  12926 

The  sample  rate  for  the  data  is  one  sample  per  sr  =  .05 ns 

The  corresponding  values  of  (d,S)  are  (.1,  .2).  With  this  choice  of  parameters,  the  forward  solve  solution 
at  2  =  0  clearly  shows  distinct  reflections  from  the  21,  22,  and  23  interfaces.  This  clear  distinction  will  aid 
in  our  approximating  the  initial  guesses,  thus  making  this  a  relatively  easy  sample  problem. 
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Figure  5:  Computed  solutions  at  different  times  of  a  windowed  electromagnetic  pulse  incident  on  a  Debye 
medium  with  a  gap. 


4.1.1  Initial  Guesses 

Assuming  the  physical  parameters  are  given  (either  known  or  from  a  previous  estimation),  we  want  to 
determine  the  depth  of  any  gap  (d),  and  the  width  of  that  gap  (d),  using  reflection  and/or  transmission 
signals.  First  we  attempt  to  get  very  close  approximations  to  d  and  6  using  information  about  the  travel 
times  of  the  data  signal,  then  we  use  these  values  as  initial  guesses  in  the  optimization  routine.  See  [3]  for  a 
complete  discussion  of  the  methods  used  in  determining  initial  estimates. 

4.1.2  Optimization  of  Least  Squares  Error 

With  initial  estimates  to  d  and  S  established,  we  define  our  inverse  problem  to  be:  find  q  =  {d,  d}  G  Qad 
such  that  the  following  least  squares  error  between  the  simulation  and  the  observed  data  is  minimized: 

1  M  S 

J(q)  =  ^5  E  EW’  #  «)  -  4)2-  (io) 

j- 1  *= i 

Here  the  Eij  are  measurements  of  the  electric  field  taken  at  M  specific  locations  (e.g.,  z°  =  0  and/or  z§  =  1) 
and  S  distinct  times  (e.g.,  every  sr  =  0.06ns).  The  E(ti,  zf;q)  are  solutions  of  the  simulations  evaluated  at 


the  same  locations  and  times  corresponding  to  the  data,  Eij,  and  using  parameter  values  q.  The  set  Qad  is 
the  feasible  set  of  q  values  determined  such  that  d  and  S  are  realistic  (e.g.  positive). 

We  apply  an  inexact  Gauss-Newton  iterative  method  to  the  optimization  problem.  That  is,  we  re-write 
the  objective  function  as 


J(q) 


1 

2  MS 


rtr 


where  Rk  =  ( E(ti ,  z°\q)  —  El;j )  for  k  =  i  +  ( j  —  1  )M  is  the  residual.  To  update  our  approximation  to  q  we 
make  the  Inexact  Newton  update  step  q+  =  qc  +  sc  where 


Sc  =  -{R\qc)TR'(qc))  *V4) 

=  -  (R'iqcfR'iqjy1  R'(qc)TR(qc) 


is  the  step,  qc  is  the  current  approximation,  and  q+  is  the  resulting  approximation.  This  is  an  inexact  method 
because  we  have  disregarded  the  S  Hessians  of  (E(ti,  z°\  q)  —  E{),  which  is  generally  acceptable  for  small 
residual  problems  [7]. 

In  this  simple  case  we  have  a  2  x  2  matrix  inverse,  so  we  can  compute  it  explicitly.  Each  iteration  requires 
one  function  evaluation  and  a  forward  difference  gradient,  which  is  two  additional  function  evaluations  (since 
we  have  two  parameters).  Each  function  evaluation  is  equivalent  to  a  simulation.  Therefore  we  want  as  few 
iterations  as  possible. 


4.1.3  Convergence 

Initial  testing  (with  2  =  0  data  only)  shows  convergence  to  8  decimal  places  of  each  parameter  in  6  iterations 
of  Gauss-Newton  with  initial  guesses  having  at  most  5%  relative  error  in  5  and  2%  in  d.  The  algorithm  does 
not  converge  to  the  correct  solution  if  the  initial  guess  for  d  has  5%  relative  error  or  if  the  initial  guess  for  6 
has  10%  relative  error. 

One  reason  that  the  algorithm  fails  to  converge  is  that  this  objective  function  is  poorly  behaved.  In 
Figure  6  we  show  a  plot  of  the  objective  function  with  respect  to  5.  The  two  very  large  peaks  in  J  on 


z2e  =  0.3, iJ  =  1 ,  N=1 024,  Nt=1 5361 ,  Ns=21 5,  MinJ(R)  =  0  @  delta  =  0.2 


delta 

Figure  6:  Nonlinear  Least  Squares  objective  function  versus  S  for  a  small  range  of  S  values  (data  at  z  =  0 
only) 
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either  side  of  the  exact  minimizer  are  due  to  the  simulated  solution  going  in  and  out  of  phase  with  the  exact 
solution.  For  this  example,  the  simulated  solution  is  most  out  of  phase  with  the  exact  solution  at  S  =  .181 
and  <5  =  .219  (i.e.,  approximately  S*  ±  |j,  which  correspond  to  the  first  and  second  peak  in  J  respectively. 
The  same  phenomenon  occurs  in  the  d  direction,  for  the  same  reasons.  See  [3]  for  a  thorough  demonstration 
of  this  behavior. 

Very  few  optimization  routines  can  provide  convergence  without  initial  conditions  between  the  two  peaks 
in  J.  The  effective  convergence  region  for  this  objective  function  applied  to  this  problem  (with  or  without 
observations  at  z  =  1)  is  within  about  8%  of  the  actual  value  of  S  when  d  is  exact  and  within  about  7.5%  of 
the  actual  value  of  d  when  5  is  exact. 

Note  also  that  the  convergence  region  is  very  dependent  on  the  frequency  of  the  interrogating  signal;  for 
higher  frequencies,  the  region  is  even  smaller.  This  is  because  the  distance  between  the  two  peaks  in  J  is 
linearly  dependent  on  A,  the  wave  length  of  the  interrogating  signal.  This  has  profound  implications  for  our 
desire  to  interrogate  with  signals  in  the  Terahertz  range.  There  are  also  peaks  in  J  with  respect  to  d  as  well, 
for  the  same  reasons. 

Further,  this  is  considering  only  a  one  parameter  minimization  problem  with  the  other  parameter  held 
fixed  at  the  exact  solution.  Convergence  is  much  worse  for  the  actual  two  parameter  problem.  In  fact,  in 
the  full  surface  plot  of  J,  a  diagonal  “trench”  occurs  approximately  along  the  line 

d=-^(S-6*)+d*  (11) 

(where  *  denotes  the  exact  solution).  This  phenomena  is  further  explained,  and  depicted  in  various  figures, 
in  [3], 

4.2  An  Improved  Objective  Function 

As  demonstrated  above,  the  usual  Nonlinear  Least  Squares  objective  function  when  plotted  with  respect  to 
either  d  or  S  has  two  large  peaks  in  J  on  either  side  of  the  exact  minimizer.  The  reason  for  these  peaks  in 
J  is  that  the  simulated  solution  goes  in  and  out  of  phase  with  the  data  as  d  or  S  change.  When  they  are 
precisely  out  of  phase,  there  is  a  very  large  absolute  error,  which  when  squared,  causes  the  objective  function 
to  have  large  peaks.  Therefore,  one  solution  is  to  not  consider  the  absolute  error,  but  instead  the  error  of 
the  absolute  values,  i.e.,  the  following  objective  function: 

.  MS  2 

Mq)  =  ^ •  (12) 

j= i *=i 

This  non-standard  mechanism  will  prevent  the  fact  that  the  signals  go  in  and  out  of  phase  with  each  other 
from  having  an  impact  on  the  objective  function,  since  positive  magnitudes  cannot  cancel  each  other  out. 
Therefore  it  gives  a  more  accurate  measure  of  the  difference  between  two  signals.  Note  that  the  orientation 
of  the  interrogating  signal  (e.g.,  peak  first)  precludes  the  possibility  of  a  solution  E  from  having  the  same 
magnitude  but  opposite  sign  as  E.  Futher,  note  that  J2(<z)  is  not  differentiable  on  a  set  of  measure  zero;  this 
is  very  unlikely  to  affect  the  finite  difference  computations  of  the  gradients,  and  did  not  present  problems  in 
our  numerical  testing.  We  plot  J2(g)  versus  5  in  Figure  7.  See  [3]  for  other  plots  including  surface  plots  as 
J2  varies  over  d  and  5. 

Note  that  while  we  have  effectively  eliminated  the  peaks  on  either  side  of  the  exact  solutions,  in  essence 
we  have  merely  converted  them  to  local  minima!  But,  since  the  minima  occur  for  the  same  reasons  the 
peaks  in  J  had  been  occurring,  they  occur  at  the  same  values  of  6.  Note  that  we  can  see  from  plotting  the 
signals  that  they  were  exactly  out  of  phase  when  they  were  shifted  by  where  A  is  the  wavelength  of  the 
interrogating  signal.  Therefore  6  is  off  by  j.  Since  we  determine  the  frequency  of  the  interrogating  signal, 
this  is  a  known  quantity,  and  we  can  predict  where  these  local  minima  will  occur  a  priori! 

Most  optimization  routines  will  continue  until  they  find  a  local  minimum,  and  since  the  two  false  minima 
described  above  are  at  least  close  to  “predictable”  locations,  we  can  easily  test  on  either  side  of  any  detected 
minimum  to  determine  if  it  is  in  fact  global.  If  J2  is  less  at  a  fixed  distance  in  the  S  direction  (e.g.,  |j  on 
either  side  of  a  detected  minimum,  i.e.,  at  either  test  location  of  the  local  minimum,  we  restart  Gauss-Newton 
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d  =  0.1,  N=1024,  Nt=1 5361 ,  Ns=215,  MinJ(R)  =  0  @  delta  =  0.2 


Figure  7:  Our  modified  Nonlinear  Least  Squares  objective  function  (J2)  versus  6  for  a  small  range  of  6  values. 
The  dotted  lines  represent  the  delta  values  that  will  be  tested  if  a  local  minimum  is  found 


at  the  new  best  guess.  Thus  if  either  of  the  two  local  minima  described  above  is  found,  we  will  eventually 
have  global  convergence  if  one  of  their  test  locations  are  sufficiently  close  to  the  global  minimizer. 

To  graphically  demonstrate  this  approach,  we  have  added  dotted  and  dash-dotted  lines  to  the  the  graph 
of  .J2(S)  in  Figure  7  to  represent  the  test  values  of  the  first  and  second  local  minima  respectively.  One  can 
clearly  see  that  if  either  local  minimum  is  detected,  the  test  value  either  to  its  right  if  it  is  the  first  one,  or 
to  its  left  if  it  is  the  second,  will  give  a  smaller  J2  and  should  eventually  lead  to  the  global  minimizer  being 
detected.  Using  this  method  we  have  in  principle  increased  our  convergence  region  to  about  25%  of  5  when 
d  is  exact.  The  same  approach  works  for  the  d  direction,  increasing  its  convergence  region  from  about  7.5% 
to  about  15%  when  S  is  exact. 

Other  possible  modifications  to  the  Least  Squares  objective  function  having  similar  effects  include  squar¬ 
ing  the  signal  instead  of  taking  the  absolute  value  (thus  preserving  smoothness  everywhere),  or  just  halving 
the  reference  signal  so  that  we  only  have  a  positive  amplitude  to  begin  with.  Each  of  these  options  will  still 
have  the  local  minima  problems  described  above,  as  well  as  their  own  unique  disadvantages. 

4.3  Testing  J2 

In  order  to  determine  the  limitations  of  an  optimization  routine  to  minimize  our  objective  function  J2  in 
a  more  practical  setting  we  examine  J2  versus  q  when  error  is  present.  In  particular  we  try  both  adding 
random  noise  to  the  data  signal,  as  well  as  testing  bad  initial  guesses  for  S  and  d.  It  should  be  noted  that 
in  the  tests  reported  on  below  we  assume  that  data  at  z  =  1  is  not  available  and  used  only  observations  at 
z  =  0. 

4.3.1  Sensitivity  to  Initial  Guesses 

For  J2  described  in  (12),  the  objective  function  is  more  sensitive  to  d  than  to  5,  therefore  it  is  imperative 
that  our  initial  guess  for  d  is  as  good  as  possible.  To  give  an  idea  of  what  may  happen  if  our  d  estimate 
were  not  within  the  15%  our  testing  has  determined  is  necessary,  we  examined  plots  of  the  objective  function 
versus  S  for  three  values  of  d,  which  are  3%,  15%,  and  30%  off  respectively  (these  are  displayed  in  [3]). 
With  errors  greater  than  15%  an  erroneous  global  minimum  appears  for  small  6  values.  This  occurs  because 
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the  first  reflection  of  the  data  is  not  matched  by  the  simulation,  but  the  second  reflection  matches  it  if  6 
is  small  enough  (see  [3]  for  details  and  sample  plots).  It  turns  out  that  the  distance  between  the  erroneous 
global  minimum  and  the  correct  minimum  is  exactly  =  0.2,  which  is  what  would  be  expected.  However, 
we  cannot  apply  the  same  idea  as  before  where  we  add  or  subtract  a  fixed  amount  to  test  for  other  local 
minima,  since  for  one,  the  “more  optimal”  of  the  two  is  farther  from  the  “true”  solution,  and  also,  we  would 
have  to  know  S  in  order  to  add  or  subtract  it  (but  S  is  what  we  are  trying  to  estimate!). 

4.3.2  Random  Observation  Noise 

In  order  to  test  the  feasibility  of  this  procedure  as  an  estimation  method,  we  have  produced  synthetic  data 
for  our  observations  £).  In  an  actual  experiment,  one  must  assume  that  the  measurements  are  not  exact.  To 
simulate  this  we  have  added  random  noise  to  the  original  signal.  The  absolute  value  of  the  noise  is  relative 
to  the  size  of  the  signal.  If  £)  is  the  data  sampled,  then  we  define  £)  =  Ei(l  +  ur]i),  where  r)i  are  independent 
normally  distributed  random  variables  with  mean  zero  and  variance  one.  The  coefficient  v  determines  the 
relative  magnitude  of  the  noise  as  a  percentage  of  the  magnitude  of  Ei,  in  particular,  v  =  0.05  corresponds 
to  10%  noise  and  v  =  0.025  to  5%  noise. 

Plots  of  the  resulting  objective  functions  for  various  values  of  v  ranging  from  2%  to  40%  are  shown  in 
[3].  Summarizing  these  results,  we  note  that  the  structure  of  the  curves  is  not  significantly  affected,  nor 
is  the  location  of  the  global  minimum.  However  the  magnitude  of  the  minimum  of  the  objective  function 
is  increased,  making  Inexact  Newton  methods  slightly  less  reliable  due  to  the  larger  residual.  Still,  our 
results  show  that  the  correct  minima  were  consistently  found  and  within  a  reasonable  amount  of  time. 
Select  examples  are  summarized  in  Table  1.  Corresponding  intial  estimates  ranged  from,  in  the  v  =  0  case, 
(d0,  S0)  =  (0.093689,0.20986)  to,  in  the  v  =  .2  case,  (d0,60)  =  (0.109668,0.172385). 


V 

d 

6 

J 

Iterations 

CPU  Time  (s) 

0 

0.1 

0.2 

1.32319E-10 

7 

160 

0.01 

0.099994 

0.199969 

0.00792792 

8 

186 

0.05 

0.099974 

0.199835 

0.199489 

13 

291 

0.2 

0.099928 

0.199204 

3.04619 

20 

435 

Table  1:  Number  of  Iterations  and  CPU  Time  for  Gauss-Newton  given  various  relative  magnitudes  of  random 
error 


5  Problem  2 

We  next  apply  the  most  useful  techniques  obtained  from  investigations  of  Problem  1  to  a  new  formulation 
of  the  interrogation  problem.  In  Problem  2  we  consider  a  dielectric  slab  and  a  metallic  backing  (conductor) 
with  a  possible  gap  between  the  two  (see  Figures  4  and  8).  Applications  of  this  specific  formulation  included 
detecting  delamination  of  insulation  from  metallic  containers,  e.g.,  insulating  foam  on  a  space  shuttle  fuel 
tank.  In  order  for  this  numerical  approach  to  be  useful  in  this  particular  application  we  must  be  able  to 
resolve  a  gap  of  width  .2 mm  inside  of  a  slab  with  a  thickness  of  at  least  20c?n  using  a  frequency  of  100 GHz. 

We  will  again  assume  the  same  physical  parameters  for  our  dielectric  and  consider  the  gap  as  a  vacuum. 
The  variables  d  and  <5  are  still  the  depth  and  the  width  of  the  gap  respectively.  One  major  difference  is  that 
in  this  problem  we  are  only  able  to  detect  the  electromagnetic  signal  in  front  of  the  material.  Also,  since  the 
metallic  backing  reflects  much  of  the  signal,  we  have  considerably  more  overlapping  of  the  reflections  to  worry 
about.  These  properties  contribute  to  the  fact  that  this  formulation  leads  to  a  much  more  difficult  inverse 
problem.  For  this  reason  we  will  be  using  more  sophisticated  optimization  routines  including  a  Levenberg- 
Marquardt  parameter  and  Implicit  Filtering.  We  will  also  need  to  develop  different  approximation  methods 
for  our  initial  guesses. 

The  implementation  of  this  problem  has  several  minor  differences  from  the  previous  one.  First,  we  now 
only  need  to  represent  two  interfaces  ij  and  Zi,  with  £q  and  £3  being  the  front  and  back  computational 
boundaries,  respectively.  Thus  now  we  define  the  depth  of  the  gap  as  d  :=  £2  —  if  and  the  width  as 
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S  :=  £3  —  £2-  Also,  as  previously  mentioned,  the  conductive  metal  backing  reflects  the  signal,  and  hence  we 
must  change  our  absorbing  boundary  conditions  at  z  =  1  (for  a  finite  computational  domain),  to  an  actual 
fixed,  Dirichlet  boundary  condition  ( E  =  0).  We  must  modify  our  finite  element  matrices  accordingly,  as 
well.  Otherwise,  the  numerical  method  for  simulation  is  the  same  as  it  was  for  Problem  1,  namely  standard 
finite  element  methods  for  spatial  derivatives,  and  an  alternating  implicit/explicit  centered  difference  time 
stepping  scheme.  Sample  solutions  are  plotted  in  Figure  9. 


Figure  8:  The  domain  of  the  material  slab  with  a  gap  between  the  medium  and  a  metallic  conductive  backing: 
f l  =  {z\zi  <  z  <  Z2}  ■ 


We  again  define  our  inverse  problem  to  be:  find  q  :=  {d,  <5}  €  Qad  such  that  an  objective  function 
representing  the  error  between  the  simulation  and  the  observed  data  is  minimized: 

min  J(q). 

ad 


Here  the  measurements  of  the  electric  field,  £),  are  taken  only  at  z  =  0,  but  still  at  S  distinct  times  (e.g., 
every  0.06ps).  The  solutions  of  the  simulations,  E(ti,0\q),  are  evaluated  at  the  same  location  and  times 
corresponding  to  the  given  data,  and  using  parameter  values  q.  In  lieu  of  actual  data  from  experiments,  we 
again  create  our  observed  data  by  using  the  simulator,  however,  the  only  information  that  is  given  to  the 
minimizer  is  the  data  observed  at  z  =  0,  which  we  will  denote  by  E. 

The  system  that  we  use  to  model  the  propagation  of  the  electric  field,  and  thus  simulate  in  order  to  solve 
our  inverse  problem,  is  as  follows,  and  includes  the  above  mentioned  Dirichlet  condition  at  z  =  1: 

Moeo(l  +  (eoo  —  1  )In)E  +  HoInP  +  —  E"  =  Js  hr  H  U 

TP  +  P  =  re0(es  -  )E  in  H 

[E  -  cE'}z~o  =  0 
[E}z= 1  =  0 
E(0,z)  =  0 
E(0,z)  =  0 

with 

Js(t,z)  =  S(z)sin(wt)I[0ttf](t). 

See  Section  2  for  a  complete  description. 
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Figure  9:  Computed  solutions  at  different  times  of  a  windowed  electromagnetic  pulse  incident  on  a  Debye 
medium  with  a  gap  between  the  medium  and  a  metallic  conductive  backing.  The  width  of  the  slab  is 
d  =  .02 m  and  the  width  of  the  gap  is  <5  =  .0002  (barely  visible  at  the  far  right  of  the  gray  region). 


5.1  Objective  Function 

As  in  the  previous  problem,  we  encounter  difficulties  when  attempting  to  use  the  standard  Least  Squares 
objective  function  to  compute  the  error  between  the  simulated  signal  and  the  observed  data.  The  constructive 
interference  of  peaks  and  troughs  produces  peaks  in  J  in  the  objective  function  on  all  sides  of  the  global 
minimum  which  make  it  nearly  impossible  to  find  the  solution  in  the  middle.  The  peaks  in  J  are  clearly 
apparent  in  Figure  10.  In  contrast,  Figure  11  shows  a  surface  plot  of  our  modified  least  squares  objective 
function 

§ 

i=  1 

It  is  clear,  as  before,  that  the  initial  guess  is  crucial  to  the  success  of  any  optimization  routine.  Notice 
that  although  J2  does  not  exhibit  the  familiar  peaks  in  J  of  Ji,  it  does  however  still  have  many  local  minima, 
which  are  just  as  difficult  to  avoid  in  a  minimization  routine. 

The  local  minima  in  J2  for  this  problem  occur  approximately  every  ^  along  the  line 


d  = 


1 


6  +  b. 
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This  happens  for  the  same  reason  as  in  Problem  1  (see  [3]  for  details  and  illustrations).  Because  we  cannot 
eliminate  these  local  minima,  we  must  appeal  to  the  procedure  that  worked  in  the  previous  problem,  namely 
testing  “check  points”.  Since  we  know  where  these  local  minima  are  occurring  with  respect  to  the  global 
minimum,  if  our  minimization  routine  finds  what  it  suspects  to  be  a  local  minima,  say  (di,<5i),  we  simply 
check  {d\  ±  <5i  =F  Q^/£ooj)i  where  a  =  1/yT  +  If  we  find  a  lower  objective  function  value,  we  restart 

our  optimization  routine  at  that  “check  point”. 


ji 


10' 


Figure  10:  Close  up  surface  plot  of  Least  Squares  objective  function  demonstrating  peaks  in  J,  and  exhibiting 
many  local  minima. 


depth 


5.2  Initial  Guesses 

In  spite  of  our  faith  in  the  “check  point”  method,  we  still  desire  to  find  the  best  initial  guesses  for  our 
optimization  routine  as  possible  so  that  we  may  hopefully  find  the  global  minimum  without  restarting.  As 
before,  we  use  the  travel  time  of  the  first  trough  to  approximate  the  location  of  the  first  interface.  However, 
in  this  formulation  we  can  take  advantage  of  some  of  the  characteristics  of  the  signals.  For  example,  the 
first  reflection  off  of  the  gap  is  always  trough-first,  and  the  second  (as  well  as  each  subsequent  reflection) 
is  always  peak-first.  For  this  reason,  if  we  want  to  locate  the  first  trough  we  can  simply  find  the  largest 
peak  (belonging  to  the  second  reflection)  and  back  track.  It  is  a  very  simple  matter  to  find  a  maximum  or 
minimum  of  a  vector  of  values.  After  the  location  of  the  largest  peak  is  found,  we  back  track  to  find  the 
minimum  in  front  of  it,  namely  that  belonging  to  the  first  reflection  off  of  the  gap.  Then  using  the  procedure 
described  in  Section  4.1,  we  approximate  the  root  immediately  in  front  of  this  trough.  That  gives  us  the 
travel  time  for  the  first  reflection  off  of  the  gap,  which  in  turn  gives  us  the  depth  d  of  the  gap. 

Finding  5  is,  unfortunately,  not  nearly  as  straightforward.  There  are  two  main  possibilities,  and  therefore, 
two  differing  approaches  to  approximating  <5,  depending  on  the  nature  of  the  reflected  signal.  We  consider 
the  two  cases: 

(i)  The  leading  trough  of  the  first  reflection  and  the  second  reflection  are  disjoint  (i.e.,  5  >  |).  In  this 
case  we  can  find  the  locations  of  the  peak  and  trough  and  use  the  travel  time  between  the  two  to 
approximate  S.  We  denote  this  approximation  by  <5i .  (Note  that  the  observed  peak  is  not  necessarily 
the  same  as  the  original  peak  unless  6  >  but  it  is  still  a  good  approximation).  See  Figure  12. 
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Figure  11:  Close  up  surface  plot  of  modified  Least  Squares  objective  function  demonstrating  lack  of  peaks 
in  J,  but  exhibiting  many  local  minima. 


(ii)  The  second  reflection  partially  truncates  the  trough  of  the  first  (i.e.,  <5  <  |).  Asa  rough  approximation, 
we  can  assume  that  the  location  of  the  actual  minimum  (trough)  is  where  the  two  signals  begin  to 
interfere  with  each  other  (the  observable  minimum).  See  Figure  13.  We  denote  this  approximation  by 

64. 

A  more  accurate  method  is  to  use  triangles  to  approximate  the  two  reflections.  By  knowing  the 
location  of  the  maximum  and  minimum  (peak  and  trough,  respectively),  and  also  the  beginning  of 
the  first  signal  (from  Section  4.1)  and  the  rough  approximation  to  the  beginning  of  the  second  signal 
using  64,  we  can  estimate  the  slopes  of  the  two  triangles  with  finite  differences.  Also  note  that  since 
the  two  signals  are  added,  the  observed  root  between  the  peak  and  trough  in  the  combined  signal  is 
actually  an  equilibrium  point  between  the  two  signals.  By  setting  equal  to  each  other  the  two  linear 
approximations  for  each  of  the  two  signals,  evaluated  at  the  equilibrium  point,  we  can  solve  for  the 
distance  between  the  starting  point  of  each  signal,  and  thus  for  S3.  See  Figure  14.  Specifically,  let 
(pi,qi)  be  the  location  of  the  trough  of  the  combined  signal  and  {p2,q2)  be  the  location  of  the  peak. 
Let  rq  be  the  location  of  the  root  in  front  of  the  trough,  and  r2  be  the  root  between  the  trough  and 
peak.  Estimate  the  slope  of  the  first  signal,  mi  <  0,  using  (pi,qi)  and  rq.  Now  if  we  let  y  =  r2  —  rq 
and  say  x  is  the  actual  distance  between  rq  and  the  beginning  of  the  second  signal,  then  setting  the 
linear  approximations  equal  in  magnitude,  but  opposite  in  sign,  at  r2  yields 

-mxy  =  m.2{y  -  x). 

Now  we  can  estimate  the  slope  of  the  second  signal,  m2  >  0,  using  {p2lq2)  and  ( r2 ,  —miy).  Also,  we 
can  re-write  the  above  equation  as 

f  —mi  +  m2\ 

x  =  -  y. 

\  m2  J 

To  find  S3  we  simply  divide  a;  by  2  and  the  (scaled)  speed  of  light  in  the  material,  i.e.,  y/e^. 

Since  each  of  the  two  situations  above  is  dependent  on  the  parameter  it  is  approximating,  we  must  also 
determine  which  of  the  above  methods  is  most  appropriate  to  use.  Thus  we  use  the  most  precise  of  the 
available  methods  to  determine  the  situation,  i.e.,  (54,  instead  of  (53  since  in  general  £3  underestimates  S  so  we 


16 


Figure  12:  The  top  plot  represents  several  signals  which  may  be  observed  in  a  simulation  of  Problem  2.  The 
bottom  plot  shows  the  sum  of  the  top  signals.  The  peak  of  the  second  signal  is  just  beginning  to  be  obscured 
by  the  first  when  <5  becomes  less  than  .  Thus  the  observable  maximum  is  still  a  good  approximation  of 
the  peak  of  the  second  signal,  and  a  trough  to  peak  distance  can  be  used  to  estimate  <5. 


do  not  want  to  use  it  as  a  criterion  for  determining  whether  <5  is  small.  (Note  that  when  <5  is  indeed  small, 
<53  is  more  accurate  than  £4.)  The  estimate  for  S 4  tends  to  be  an  overestimate,  and  is  only  valid  if  S  <  . 

Unfortunately,  (5i  also  tends  to  be  an  overestimate,  so  we  prefer  to  only  trust  it  entirely  if  it  is  larger  than 
j.  If  neither  (5i  nor  (>3  is  a  sufficient  approximation  we  choose  to  use  the  average  of  the  two,  and  call  it  82- 

Therefore  our  algorithm  for  approximating  <5  is  as  follows: 

(a)  If  ^4  <  -|  then  use  (53 

(b)  else  if  (5i  >  |  then  use  (5i 

(c)  else  use  82  (average  between  (5i  and  S3). 

We  tested  our  approximating  methods  on  exact  depth  ( d )  values  of:  .02,  .04,  .08,  .1,  and  .2  to,  and 
values  of  width  (5):  .0001,  .0002,  .0004,  .0006,  and  .0008  to.  Since  |  is  the  transition  point  between  the 
two  situations,  it  is  understandable  why  S  close  to  this  value  is  the  most  difficult  to  accurately  resolve.  We 
chose  this  range  of  S’ s  because  our  choice  of  frequency  gives  |  =  3.7475  x  10_4m.  See  the  tables  in  [3]  for 
the  initial  estimates  of  d  and  5. 

The  approximations  improve  slightly  as  the  number  of  finite  elements  is  increased,  and  appeared  to 
converge  to  fixed  values.  This  suggests  that  numerical  error  (and  instability)  can  affect  the  estimates.  For 
each  case  there  is  a  significant  amount  of  visible  numerical  error  in  the  simulations  below  a  certain  number 
of  elements,  therefore  in  approximating  S  we  chose  to  use  the  number  of  elements  just  above  the  threshold. 

While  the  initial  estimates  were  relatively  inaccurate,  some  8  approximations  being  almost  100%  off  from 
the  true  solution  value,  in  the  numerical  tests  we  performed,  all  the  initial  estimates  were  sufficiently  close 
to  the  true,  global  minima  as  to  not  cause  the  optimization  routine  to  result  in  a  false,  local  minima.  While 
our  “check  point”  method  is  available  if  needed,  it  is  much  more  efficient  to  have  an  accurate  initial  estimate 
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Figure  13:  The  top  plot  represents  several  signals  which  may  be  observed  in  a  simulation  of  Problem  2.  The 
bottom  plot  shows  the  sum  of  the  top  signals.  The  trough  of  the  first  signal  is  partially  truncated  by  the 
second  signal.  In  this  case  the  observed  minimum  is  a  still  a  good  approximation  to  where  the  second  signal 
begins.  For  smaller  S,  a  linear  approximation  must  be  used. 


than  to  restart  after  optimizing  from  a  bad  one.  Still,  the  report  [3]  describes  several  very  real  examples 
where  the  “check  point”  method  would  be  a  necessary  last  resort. 

5.3  Optimization  Method 

Now  that  we  have  approximated  our  initial  guesses,  we  need  to  minimize  the  objective  function  in  order  to 
solve  the  inverse  problem.  In  Problem  1,  Gauss-Newton  was  sufficient  to  find  the  global  minimum  for  most 
cases.  In  this  formulation,  however,  we  will  apply  more  sophisticated  methods,  reverting  to  Gauss-Newton 
whenever  possible  since  its  convergence  rate  is  best. 

The  first  modification  we  make  to  Gauss-Newton  is  to  add  a  Levenberg-Marquardt  parameter,  vc  (see 
[7]).  The  Inexact  Newton  step  becomes 

sc  =  -  (R1  (qc)T R' (qc)  +  ucI )  1  R'(qc)TR(qc). 

The  parameter  adds  regularization  by  making  the  model  Hessian  positive  definite.  The  method  uses  a 
quadratic  model  Hessian,  and  also  has  a  built-in  line  search  with  a  sufficient  decrease  condition.  The  line 
search  is  based  on  the  predicted  decrease  computed  from  the  quadratic  model.  If  the  actual  improvement 
of  the  objective  function,  J,  is  close  to  the  amount  predicted  by  the  model  Hessian  after  a  step  is  taken, 
then  the  method  decreases  the  Levenberg-Marquardt  parameter,  uCl  effectively  increasing  the  relative  size 
of  the  next  step,  which  hopefully  accelerates  the  convergence.  As  vc  is  decreased  to  0  the  method  becomes 
Damped  Gauss-Newton  (meaning  Gauss-Newton  with  a  line  search).  If,  however,  the  actual  improvement 
of  J  after  a  step  is  not  sufficient  (or  is  even  negative),  vc  is  increased,  effectively  scaling  back  the  Newton 
step,  and  we  retest.  If  there  are  too  many  reductions  then  we  declare  a  “line  search  failure”  meaning  that 
too  small  a  step  is  required  to  decrease  the  objective  function. 
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Figure  14:  This  schematic  shows  the  roots,  extrema,  distances,  and  slopes  used  in  the  computation  of  £3. 


Usually  a  method  would  exit  after  a  line  search  failure,  returning  the  best  approximation  so  far.  But 
we  use  this  failure  to  call  an  adaptive  mesh  size  routine,  i.e. ,  an  Implicit  Filtering  technique.  The  idea  is 
that  the  failure  is  likely  due  to  the  fact  that  the  direction  the  finite  difference  gradient  chose  is  probably 
not  an  actual  “descent  direction”  in  the  global  sense.  In  other  words,  the  finite  differencing  is  most  likely 
differentiating  noise.  In  the  same  manner  that  a  smooth  surface  may  look  rough  under  a  microscope,  using 
too  small  of  a  differencing  step  amplifies  effects  from  round-off  error  and  other  sources  of  numerical  noise. 
Our  technique  is  to  increase  the  relative  differencing  step,  h,  recompute  the  gradients,  and  then  try  the 
Levenberg-Marquardt  method  again.  The  relative  differencing  step,  h,  is  such  that  the  gradient,  V^,  of 
J(q)  =  J([d,  <J])  is  computed  with 


V^([d,<S]) 


hd 

j(d,(l+h)6)-j(d,6) 

hS 


We  apply  a  similar  approach  to  modifying  the  differencing  step  h  as  we  do  for  changing  vc  in  that  after  a 
successful  step  we  decrease  h ,  but  if  we  have  another  failure  we  increase  h  even  more.  Since  the  convergence 
rates  of  gradient  based  methods  are  dependent  on  the  size  of  h  (for  example  Gauss-Newton  is  0(h2)),  we 
want  h  to  be  as  small  as  possible  and  still  be  effective,  similarly  with  vc.  We  use  a  three  tiered  approach 
to  changing  h.  Initially  we  set  h  =  1CU9.  To  increase  h  we  raise  it  to  the  |  power,  to  decrease  we  raise  it 
to  the  |  power.  Additionally  we  define  1CU4  to  be  the  maximum  allowable  differencing  step  value.  Thus 
h  G  {10-9,10_6,10-4}. 

In  general  an  optimization  method  exits  with  “success”  if  the  norm  of  the  current  gradient  is  less  than 
tol  times  the  norm  of  the  initial  gradient.  However,  in  our  method  we  do  not  immediately  trust  the  finite 
difference  gradients,  and  instead  call  Implicit  Filtering  again  when  the  gradients  appear  small.  When  we 
have  verified  small  gradients  on  all  three  scales  (the  various  values  of  the  differencing  step  h  defined  above), 
then  we  exit  with  “success” . 
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Remark  2  In  practice,  a  very  good  solution  is  found  within  a  couple  of  Levenberg-Marquardt  steps,  and 
then  an  equal  number  of  Implicit  Filtering  iterations  verify,  and  sometimes  enhance,  this  solution.  In  the 
interest  of  efficiency,  and  since  this  is  a  parameter  identification  problem,  we  exit  early  with  “success”  if  our 
objective  function  is  satisfactorily  small  (i.e.,  tol  times  the  initial  value),  which  can  save  on  average  about 
half  of  the  possible  iterations. 

Additionally  we  impose  a  restriction  on  the  number  of  “pullbacks”  on  each  linesearch,  and  on  the  number 
of  iterations,  effectively  limiting  the  total  number  of  function  calls.  If  a  small  gradient  has  not  been  verified 
on  all  scales  before  exhausting  the  maximum  number  of  iterations,  we  exit  with  “failure” . 

5.4  Numerical  Issues 

For  small  N  the  difficult  cases  are  those  with  large  depth.  This  is  because  the  computational  domain  is 
effectively  increased  when  the  depth  is  increased,  making  the  mesh  sizes  larger  and  increasing  the  level  of 
numerical  error.  The  magnitude  of  S  does  not  seem  to  have  a  significant  effect  on  the  convergence  of  the 
method. 

An  obvious  disadvantage  to  having  a  large  N  is  that  each  simulation  takes  much  longer.  In  general  the 
total  execution  time  is  quadrupled  when  the  number  of  elements  is  doubled.  This  is  consistent  with  the 
fact  that  complexity  of  the  most  time  consuming  part  of  the  simulation,  the  linear  solves,  is  O(N),  and  the 
number  of  time  steps  Nt  is  also  O(N).  So  when  we  double  the  number  of  finite  elements  we  are  also  doubling 
the  number  of  time  steps.  Therefore,  we  get  an  overall  complexity  of  0(N2).  Thus,  as  mentioned  before,  in 
our  inverse  problem  we  choose  to  use  the  number  of  elements  just  above  the  threshold  of  when  numerical 
error  is  apparent. 

We  should  also  mention  that  in  order  to  create  data,  in  lieu  of  actual  experimental  data,  we  perform  a 
simulation  at  a  higher  resolution  believing  it  to  be  more  accurate.  Specifically,  we  double  the  number  of 
finite  elements.  Since  the  time  step,  and  therefore  the  effective  sample  rate  if  the  time  step  is  too  large,  are 
both  dependent  upon  the  mesh  size  (refer  to  [3]),  the  sample  times  of  the  simulated  data  do  not  necessarily 
correspond  with  the  sample  times  of  the  simulations  at  the  lower  resolution.  (In  general  we  have  twice 
as  many  samples  from  the  higher  resolution.)  Thus  in  order  to  compute  the  modified  least  squares  error 
between  the  two  vectors,  we  perform  a  linear  interpolation  of  the  simulated  data  onto  the  sample  times  at 
the  lower  resolution.  See  Figure  15.  Note  that  in  the  usual  case  where  we  simply  have  twice  as  many  sample 
points  from  the  higher  resolution  simulation,  we  are  in  effect  discarding  sample  points  rather  than  doing  a 
true  interpolation. 

For  comparison  we  compute  the  low  resolution  simulation  using  the  values  d*  and  d*  (note  that  this  is 
not  the  same  as  taking  the  high  resolution  simulation  and  interpolating  it  onto  the  low  resolution  time  steps, 
which  we  actually  use  as  our  observed  data).  In  every  case  that  we  have  tested,  J,  when  computed  with 
the  d  and  6  values  found  from  the  optimization  routine  ( dmin  and  Smln),  is  less  than  or  equal  to  J  when 
computed  with  the  original  values  (d*  and  5*).  This  suggests  that  an  actual  global  minimum  of  the  objective 
function  has  been  found,  even  though  the  final  estimates  of  d  and  5  themselves  are  not  necessarily  equal 
to  d*  and  d*.  Note  in  Figure  15  that  the  simulation  using  original  values,  ( d*,6 *),  is  in  fact  closer  to  the 
original  data,  but  the  simulation  using  the  minimizer  values,  (dmin,  Smin),  is  closer  to  the  interpolated  data 
(see  for  example  the  [.335,  .3352]  interval). 

Although  we  could  compute  our  optimization  routine  at  the  same  resolution  as  the  simulated  data  to 
get  a  better  fit  in  our  tests,  this  would  not  properly  represent  the  real-life  phenomenon  of  sampling  data. 
Sampled  data  is  inherently  not  a  completely  accurate  representation  of  a  physical  observation.  We  believe 
that  our  interpolation  approach  gives  a  more  realistic  expectation  of  how  our  method  would  perform  given 
actual  experimental  data.  In  order  to  further  test  the  robustness  of  our  inverse  problem  solution  method  we 
introduce  random  noise  to  the  detected  data  in  Section  5.5. 

5.5  Numerical  Results 

Tables  2  and  3  show  the  final  computed  approximations  for  the  depth  of  the  slab  (drnm  )  and  the  width  of 
the  gap  behind  it  (5min).  The  relative  differences  from  the  original  values  used  to  generate  the  data  (d*  and 
5*),  are:  for  depth,  on  the  order  of  .0001  and  for  5,  on  the  order  of  .01.  However,  this  does  not  imply  that 
the  optimization  routine  was  unable  to  find  the  optimal  solution.  Recall  that  since  our  data  is  generated 
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Figure  15:  Plotted  are  the  actual  simulated  data  ( N  =  2048),  the  interpolation  of  the  simulated  data  onto 
the  low  resolution  sample  times  (N  =  1024),  the  result  of  the  minimization  routine  (N  =  1024),  and  a  low 
resolution  (N  =  1024)  simulation  using  the  exact  values  of  d  and  5. 


with  essentially  a  different  simulator  than  our  forward  solves,  the  original  values  do  not  necessarily  minimize 
the  objective  function.  The  objective  function  values  give  a  better  indicator  of  how  well  the  optimization 
routine  works  since  it  shows  the  fit  to  the  generated  data.  Table  4  shows  the  final  objective  function  values. 
In  each  of  these  cases,  the  final  objective  function  value  (, Jmin )  was  less  than  J*  :=  J(q*).  In  fact,  the  ratios 
Jr  :=  Jmin/ J*  were  on  average  .3008.  We  consider  any  Jr  <  1  to  represent  a  successful  convergence. 

Although  6  values  that  are  near  |  =  3.7475  x  10_4m  are  the  most  difficult  for  which  to  obtain  initial 
approximations,  we  see  that  the  objective  function  values  in  these  cases  are  just  as  small  (and  the  final 
estimates  are  just  as  close)  as  for  other  5  values. 

The  execution  time,  in  seconds,  as  well  as  the  number  of  function  calls,  are  given  in  [3].  While  the  above 
tables  establish  that  we  were  actually  able  to  resolve  the  case  of  20cm  depth,  there  was  a  price  we  had  to  pay. 
The  average  execution  times  for  each  of  different  mesh  sizes  ( N  =  1024,2048,4096,8192,  and  16384)  were 
39,  248, 1452,  6229,  and  35509  seconds,  respectively.  Each  represents  an  increase  in  time  over  the  previous 
mesh  size  by  a  factor  of  6.4,5.9,4.28,  and  5.7,  respectively.  This  is  consistent  with  the  fact  that  the  forward 
solves  are  order  0(h2).  However,  the  additional  sample  points  for  the  larger  N  cases  allowed  for  smaller 
initial  objective  function  values  which  resulted  in  increasingly  more  iterations  to  satisfy  the  relative  tolerance 
in  our  stopping  criteria.  This  explains  why  we  do  not  see  ratios  closer  to  the  expected  4  for  order  0(h2) 
methods. 

5.5.1  Relative  Random  Noise 

We  add  random  noise  to  the  signal,  as  mentioned  above,  in  order  to  more  closely  simulate  the  experimental 
process  in  data  collection.  As  in  Section  4.3.2,  we  start  with  relative  noise  where  the  absolute  value  of  the 
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6 


d 

.0001 

.0002 

.0004 

.0006 

.0008 

.02 

(N=1024) 

0.0200053 

0.0200022 

0.0200006 

0.0200005 

0.0200002 

.04 

(N=2048) 

0.0399948 

0.0399974 

0.0400005 

0.0400005 

0.0399999 

.08 

(N=4096) 

0.0799973 

0.0799987 

0.0800006 

0.0800006 

0.0800003 

.1 

(N=8192) 

0.0999945 

0.0999974 

0.1 

0.1 

0.0999999 

.2 

(N=16384) 

0.200011 

0.200005 

0.2 

0.2 

0.200001 

Table  2:  The  final  estimates  of  d. 


6 


d 

.0001 

.0002 

.0004 

.0006 

.0008 

.02 

(N=1024) 

9.40622e-05 

0.000196754 

0.000398642 

0.000597275 

0.00079707 

.04 

(N=2048) 

0.000106435 

0.000203916 

0.000394204 

0.000592156 

0.000793622 

.08 

(N=4096) 

0.000103585 

0.000202273 

0.000395791 

0.000593861 

0.000794401 

.1 

(N=8192) 

0.000106593 

0.000203876 

0.000396203 

0.000594976 

0.000795985 

.2 

(N=16384) 

8.7456e-05 

0.000191808 

0.00040297 

0.000602902 

0.00080129 

Table  3:  The  final  estimates  of  6. 


6 


d 

.0001 

.0002 

.0004 

.0006 

.0008 

.02 

(N=1024) 

0.00786171 

0.00906699 

0.0115657 

0.0233783 

0.0447687 

.04 

(N=2048) 

0.021516 

0.0343314 

0.0514108 

0.0700747 

0.0927117 

.08 

(N=4096) 

0.0116105 

0.0145428 

0.0201004 

0.0272513 

0.0344458 

.1 

(N=8192) 

0.00304723 

0.00547532 

0.00779186 

0.00931778 

0.0118529 

.2 

(N=16384) 

0.000609258 

0.00133978 

0.00146975 

0.000962975 

0.000766141 

Table  4:  The  objective  function  value  of  the  final  estimates. 


noise  is  proportional  to  the  size  of  the  signal.  If  Ei  is  the  data  sampled,  then  we  define  Et  =  77,(1  +  vrr}i), 
where  %  are  independent  normally  distributed  random  variables  with  mean  zero  and  variance  one.  Again, 
the  coefficient  vr  determines  the  relative  magnitude  of  the  noise  as  a  percentage  of  the  magnitude  of  Ei,  in 
particular,  vr  =  0.01  corresponds  to  2%  noise.  We  tested  relative  magnitude  levels  of  2%,  10%,  and  20% 
(corresponding  to  vr  —  .01,. 05,  and  .1  respectively).  See  [3]  for  tables  of  initial  estimates.  In  nearly  all 
the  cases  the  estimate  was  close  enough  for  the  optimization  method  to  converge  ( Jr  <  1)  to  the  expected 
minimum.  The  only  exceptions  were  with  v  =  .1  and  6  =  .0004,  which  are  understandably  the  most  difficult 
cases. 

The  final  approximations  dmin  and  Smin  in  the  presence  of  noise  are  also  given  in  [3] .  Some  approximations 
with  high  noise  appear  to  be  better  approximations  than  some  with  little  or  no  noise.  For  example,  with 
6*  =  .0001,  d*  =  .04,  the  vr  =  .1  final  approximations  are  an  order  of  magnitude  closer  to  the  original  values 
than  the  vr  =  0  final  approximations.  This  is  not  to  say  that  the  noise  helps  the  approximation  method. 
Rather,  it  is  for  the  same  reason  that,  for  example,  as  shown  in  Figure  15,  the  actual  parameter  values 
produced  a  signal  farther  away  (in  the  Least  Squares  sense)  from  the  generated  data  than  a  signal  computed 
with  the  approximated  parameter  values.  The  resulting  objective  function  values  give  a  better  indication  of 
the  accuracy  of  the  approximation  to  the  data.  The  final  objective  function  values  corresponding  to  ur  =  0 
were  two  orders  of  magnitude  smaller  on  average  than  those  resulting  from  vr  =  .1.  Thus,  it  is  clear  that 
the  data  without  noise  is  more  accurately  matched  by  its  approximations  than  those  with  noise. 

5.5.2  Standard  Error  Analysis 

In  an  actual  inverse  problem  using  data  collected  by  experiment,  one  desires  to  have  confidence  intervals 
on  all  parameter  estimates.  We  will  apply  standard  error  techniques  to  an  Ordinary  Least  Squares  (OLS) 
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formulation  of  our  problem  to  obtain  confidence  intervals  on  our  estimates.  In  order  to  rewrite  our  objective 
function  in  an  OLS  formulation,  we  define  y(t;q)  =  \E(t,0;q)\  to  be  our  estimate  to  y—  |2?|,  which  is  the 
data  we  are  trying  to  fit  by  determining  q  =  ( d1  (5).  Now  it  is  clear  that  our  objective  function  can  be  written 
in  the  standard  OLS  form 

1  Na 

J(q )  =  -yi\2  ■ 

s  i—  1 

For  simplicity  of  terminology,  in  this  section  alone,  we  will  refer  to  \Ei\  as  the  data  and  to  | E(ti,  0;  q)  as  the 
simulations. 

With  the  relative  random  noise  described  above  we  do  not  have  constant  variance,  as  is  demonstrated  in 
Figures  16  and  17.  Here  we  have  plotted  the  residual  r,  :=  \E(ti,0;qoLs)\  ~  \Ei\  against  time,  tj,  and  also 
against  \E(U,  0;  qoLs)\-  As  one  would  expect  with  noise  that  is  relative  in  size  to  the  signal  value,  we  have 
a  pattern  in  Figure  16  that  follows  the  pattern  of  the  original  signal.  Figure  17  demonstrates  the  fan  shape 
associated  with  noise  that  is  dependent  upon  the  size  of  the  signal,  i.e.,  non  constant  variance. 

Since  constant  variance  is  most  conveniently  assumed  in  standard  error  analysis,  we  further  consider 
estimates  obtained  from  an  inverse  problem  applied  to  data  with  constant  variance  random  noise  added.  In 
particular,  the  data  we  now  consider  is  generated  by 

Ej  =  E(tj,  0,  ;q*)+  f3vrVj 


where 

Vj  ~  A/"(0, 1) 

and  the  constant  (3  is  a  scaling  factor  chosen  simply  so  that  the  noise  level,  vr,  will  somewhat  correspond  to 
the  parameter  ry  used  in  the  previous  section  on  relative  noise.  Specifically,  /3  =  max;  Ei/ 10  ensures  that  J* 
in  the  constant  variance  cases  is  on  the  same  order  of  magnitude  as  those  in  the  relative  noise  cases  above 
for  all  choices  of  d  and  5  that  we  have  considered. 

The  variance  of  this  data  is 

^  =  £[p2v2rV2}  =  (32v2r£  [„?]  =  / 32u2 

where  £  denotes  the  expectation.  Therefore,  we  do  have  constant  variance.  Note  further  the  resulting  lack  of 
patterns  in  Figures  18  and  19.  The  suspicious  looking  phenomenon  of  many  points  on  the  line  E  =  0  is  simply 
because  in  the  original  data  E  is  very  close  to  zero  most  of  the  time.  Figure  20  demonstrates  graphically 
the  difference  between  relative  noise  and  constant  variance  noise.  The  relative  noise  case  is  particularly 
difficult  in  our  inverse  problem  since  most  of  our  initial  estimates  are  based  on  accurately  determining  the 
peak  locations,  yet  this  is  exactly  where  most  of  the  relative  noise  is  concentrated. 

With  constant  variance,  and  further,  assuming  that  each  ly  is  identically  independently  (normally)  dis¬ 
tributed,  we  have  that  (see  [4])  in  the  limit  as  Ns  — »  oo 

qOLS  ~  -A/2  (<7o,  ctq  [ST(q0)S(q0)]  *)  • 

Here  S(q)  =  which  is  an  Ns  x  2  matrix  since  q  =  ( d,S )  and  \E\  is  evaluated  at  Ns  sample  times. 

Also,  the  scale  parameter  is  approximately  given  by 


In  the  above  equations,  <70  denotes  the  theoretical  “true”  value  of  the  parameter  that  best  describes  the 
system  from  which  the  data  is  taken.  Note  that  in  this  case,  this  is  not  necessarily  the  same  as  q*  since  the 
method  used  to  generate  the  data  is  different  from  the  forward  solve  simulator.  Therefore  q$  is  generally 
unknown  even  in  examples  with  simulated  data. 

As  demonstrated  in  the  previous  sections,  our  qoLS  is  often  a  better  minimizer  than  even  the  original  value 
of  q*,  therefore  we  will  approximate  qo  in  the  above  equations  by  qoLS-  In  particular,  if  we  denote  the  covari¬ 
ance  matrix  as  Co  =  Uq  [St (<7o)<S(<7o)]  >  then  we  will  approximate  Co  by  C  —  o2qLS  [ST (qoLs)S{qc>Ls)]  , 

where 

1  ns  2 

aOLS  =  AT  2  ^  (\E(ti>0',qOLs)\  —  • 

s  i=  1 


23 


We  compute  a qLS  by  multiplying  our  Jmin  by  an  appropriate  conversion  factor,  since  they  are  defined  in  a 
similar  manner.  However,  in  order  to  compute  the  partial  derivatives  with  respect  to  d  and  8  in  S  we  employ 
forward  differencing,  which  requires  an  additional  forward  simulation  for  each  qj.  For  q  =  qoLS  we  have,  for 
example 

c  d\E\  \E(ti,0;[q1,q2})\-\E(ti10;[(l- hd)q1,q2])\ 

«->il  —  - \ — x - 

oqi  hdqi 

and  similarly  for  each  <Si2.  In  our  computations  we  used  the  relative  differencing  factor  of  hd  =  1  x  10-4. 
One  could  also  use  a  sensitivity  equations  approach  (e.g.,  see  [1]  and  the  references  therein),  but  since  the 
variational  equations  are  quite  difficult  to  solve  for  this  example,  we  choose  instead  to  approximate  the 
partials  with  respect  to  q  directly  with  our  simulations. 

We  also  need  to  point  out  that  while  taking  the  absolute  value  of  a  function  limits  differentiability  at  a 
small  number  of  points,  the  derivative  does  exist  almost  everywhere.  The  absolute  value  function  does  not 
change  the  magnitude  of  the  derivative  where  it  exists,  which  is  what  we  need  to  compute  the  dot  product 
of  S  with  itself.  By  using  finite  differences  to  estimate  derivatives,  we  are  essentially  under-estimating  at 
the  discontinuities.  Under-estimating  a  few  points  out  of  thousands  is  not  going  to  significantly  change  our 
covariance  matrix.  (Alternatively,  one  could  have  defined  the  objective  function  by  squaring  the  signals 
instead  of  taking  absolute  values  to  avoid  this  problem.  In  this  research  we  were  interested  in  comparing  Ji 
and  J2  in  previous  sections  above  and  changing  the  scale  of  E  by  squaring  it  would  have  prevented  this.) 

With  S  calculated,  we  can  now  evaluate  C  =  &ols  [ST  (qoLs)S(qoLs)]  ■  Then  the  standard  error  for 
qi  =  d  is  estimated  by  \JC\\  while  the  standard  error  for  q2  =  8  is  estimated  by  yJC22.  See  Tables  5  through 
12  for  confidence  intervals  relating  to  various  d*,  8*  and  vr  values.  For  example,  in  the  case  of  d*  =  .02, 
8*  =  .0002  and  with  vr  =  .01  our  covariance  matrix  is 


'  2.37122  x  10"15  -4.43815  x  10"15' 

-4.43815  x  10"15  9.1829  x  10"15 


which  results  in  the  confidence  intervals  d  €  (2.00004±4.86952  x  10  6)  x  10  2  and  8  €  (1.9941±0. 000958274)  x 
10"4. 

The  width  of  these  bounds  are  ±0.000243471%  and  ±0.0480555%  of  the  approximation  value  respectively. 
For  the  d*  =  .02  case,  the  average  size  of  the  confidence  intervals  for  vr  =  .01,  .05,  .1  respectively  were 
±.0002%,  ±.001%,  ±.002%  (averaged  over  various  <5*  values  ranging  from  .0001  to  .0008).  It  is  interesting  that 
the  widths  of  the  confidence  intervals  nearly  exactly  double,  on  average,  when  the  noise  level  doubles.  For  the 
d*  =  .04  case  the  average  size  of  the  confidence  intervals  were  ±.0001%,  ±.0006%,  ±.001%.  Likewise,  when 
the  widths  of  the  confidence  intervals  for  8*  =  .0002  are  averaged  over  several  various  d*  values  (.02,  .04,  .08) 
we  get  ±.05999%,  ±.2883%,  ±.5718%  for  vr  =  .01,  .05,  .1  respectively.  For  <5*  =  .0004  the  averages  are 
±.03331%,  ±.1575%,  ±.3154%.  In  general,  larger  d*  and  8*  values  have  smaller  (tighter)  confidence  intervals. 
This  suggests  that  the  approximations  found  in  these  cases  are  better  than  those  estimating  small  parameters. 
While  this  is  intuitive,  it  is  not  apparent  looking  at  the  estimates  themselves  or  even  the  final  objective 
function  values  (see,  for  example,  Table  4). 
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|r(q)|vs  t,  (r(q)=  |e(q)|-|ehat|) 


Figure  16:  Plots  of  the  absolute  value  of  the  residual  r;  =  \E{U,  0;  qoLs) \  —  \Ei\  versus  time  t,  when  the  data 
contains  relative  random  noise. 

|r(q)l  vs  |e(q)i,  (r(q)=  |e(q)|-|ehat|) 
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Figure  17:  Plots  of  the  absolute  value  of  the  residual  r,;  =  \E(ti,  0;  qoLs) I  —  \Ei\  versus  the  absolute  value  of 
the  electric  field  E(ti,  0;  qoLs)  when  the  data  contains  relative  random  noise. 
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Ir(q)l  vs  t,  (r(q)=  |e(q)|-|ehat|) 


t  (ns) 


Figure  18:  Plots  of  the  absolute  value  of  the  residual  =  \E(ti,  0;  qoLs) I  —  \Ei\  versus  time  ti  when  the  data 
contains  constant  variance  random  noise. 


Ir(q)l  vs  |e(q)|,  (r(q)=  |e(q)|-|ehat|) 


Figure  19:  Plots  of  the  absolute  value  of  the  residual  r*  =  | E(ti,  0;  qoLs)  I  —  \Ei\  versus  the  absolute  value  of 
the  electric  field  E(ti,0;  qoLs)  when  the  data  contains  constant  variance  random  noise. 
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electric  field  (volts/meter) 


Relative  Noise  vs  Constant  Variance 


Figure  20:  The  difference  between  data  with  relative  noise  added  and  data  with  constant  variance  noise 
added  is  clearly  evident  when  E  is  close  to  zero  or  very  large. 
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6 

.0002 

.0004 

.0008 


d*  =  .02  (N  =  2048) 

(2.00005  ±  0.30284  x  10~7)  x  10 
(2.00001  ±  6.50411  x  10"7)  x  10 
(2.00001  ±4.91232  x  10"7)  x  10 


-2 

-2 


d*  =  .04  (N  =  4096) 

(4.00013  ±  1.62162  x  10~6)  x  KT2 
(4.00001  ±  1.19064  x  10"6)  x  10"2 
(4.00002  ±  9.05240  x  10"7)  x  10"2 


Table  5:  Confidence  intervals  for  the  OLS  estimate  of  d  when  the  data  is  generated  with  no  noise  (i.e. 
vr  =  0.0). 


<5 

.0002 

.0004 

.0008 


d*  =  .02  (N  =  2048) 

(2.00004  ±4.86952  x  10~B)  x  10 
(2.00001  ±  3.50259  x  10"6)  x  10 
(2.00001  ±  2.87772  x  10"6)  x  10 


-2 

-2 


d*  =  .04  (N  =  4096) 

(4.00013  ±  5.69385  x  10~H)  x  KT2 
(4.00001  ±  4.02428  x  10"6)  x  10"2 
(4.00001  ±  3.32933  x  10"6)  x  10"2 


Table  6:  Confidence  intervals  for  the  OLS  estimate  of  d  when  the  data  is  generated  with  noise  level  vr  =  .01 


6 

.0002 

.0004 

.0008 


d*  =  .02  (N  =  2048) 

(2.00U04  ±2.41541  x  10“5)  x  10 
(2.00000  ±  1.68896  x  KT5)  x  10 
(2.00003  ±  1.40398  x  10"5)  x  10 


-2 

-2 


d*  =  .04  (N  =  4096) 

(4.00014  ±  2.76640  x  10~5)  x  10”2 
(4.00001  ±  1.90853  x  10"5)  x  10"2 
(4.00000  ±  1.60390  x  10"5)  x  KT2 


Table  7:  Confidence  intervals  for  the  OLS  estimate  of  d  when  the  data  is  generated  with  noise  level  vr  =  .05 


6 

.0002 

.0004 

.0008 


d*  =  .02  (N  =  2048) 

(2.00000  ±4.72903  x  KT5)  x  10 
(2.00003  ±  3.39327  x  KT5)  x  10 
(2.00003  ±2.79911  x  10"5)  x  10 


-2 

-2 


d*  =  .04  (N  =  4096) 

(4.00014  ±  5.48283  x  10~5)  x  10"2 
(4.00002  ±  3.87474  x  10"5)  x  10"2 
(4.00003  ±  3.19526  x  10"5)  x  10"2 


Table  8:  Confidence  intervals  for  the  OLS  estimate  of  d  when  the  data  is  generated  with  noise  level  ur  =  .  1 
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s 

.0002 

.0004 

.0008 


d*  =  .02  (N  =  2048) 

(1.99272  ±  0.000182978)  x  10 
(4.00035  ±  0.000201885)  x  10 
(7.99833  ±  0.000136586)  x  10 


-4 

-4 


d*  =  .04  (N  =  4096) 

(1.98142  ±  0.000317616)  x  10~4 
(4.00737  ±  0.000369841)  x  10"4 
(8.00332  ±  0.000251291)  x  10"4 


Table  9:  Confidence  intervals  for  the  OLS  estimate  of  6  when  the  data  is  generated  with  no  noise  (i.e. 
vr  =  0.0). 


S 

.0002 

.0004 

.0008 


d*  =  .02  (N  =  2048) 

(1.99410  ±  0.000958274)  x  10~4 
(4.00170  ±  0.00108740)  x  10"4 
(7.99882  ±  0.000800042)  x  10“4 


d*  =  .04  (N  =  4096) 

(1.98029  ±  0.00111475)  x  10~4 
(4.00667  ±  0.0012499)  x  10"4 
(8.00486  ±  0.000923838)  x  10"4 


Table  10:  Confidence  intervals  for  the  OLS  estimate  of  <5  when  the  data  is  generated  with  noise  level  vr  =  .01 


S 

.0002 

.0004 

.0008 


d*  =  .02  (N  =  2048) 
(1.99606  ±  0.00475672)  x  10 
(4.00190  ±  0.00524360)  x  10 
(7.99045  ±  0.00391181)  x  10 


-4 

-4 


d*  =  .04  (N  =  4096) 

(1.98106  ±  0.00541764)  x  10~4 
(4.01214  ±  0.00593246)  x  10"4 
(8.00947  ±  0.00444525)  x  10“4 


Table  11:  Confidence  intervals  for  the  OLS  estimate  of  6  when  the  data  is  generated  with  noise  level  vr  =  .05 


5 

.0002 

.0004 

.0008 


d*  =  .02  (N  =  2048) 

(2.00017  ±  0.00932701)  x  10~4 
(4.00070  ±  0.0105331)  x  10"4 
(7.99698  ±  0.00778563)  x  10"4 


d*  =  .04  (N  =  4096) 

(1.97674  ±  0.0107203)  x  10~4 
(4.01229  ±  0.0120445)  x  10"4 
(8.00361  ±  0.00886925)  x  10"4 


Table  12:  Confidence  intervals  for  the  OLS  estimate  of  <5  when  the  data  is  generated  with  noise  level  vr  =  .1 
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6  Conclusion 


In  this  presentation,  we  have  explored  a  “proof  of  concept”  formulation  of  an  inverse  problem  to  detect  and 
characterize  voids  or  gaps  inside  of,  or  behind,  a  dielectric  medium.  We  have  simplified  the  problem  to  one 
dimension  and  used  Maxwell’s  equations  to  model  a  pulsed,  normally  incident  electromagnetic  interrogating 
signal.  We  use  Finite  Element  discretization  in  space,  and  Finite  Differences  in  time,  to  simulate  the  electric 
field  in  the  time  domain.  This  is  coupled  with  a  Levenberg-Marquardt  scheme  in  an  optimization  step  with  an 
innovative  cost  functional  appropriate  for  reflected  waves  where  phase  differences  can  produce  ill-posedness 
in  the  inverse  problem  when  one  uses  the  usual  ordinary  least  squares  criterion.  We  have  successfully 
demonstrated  that  it  is  possible  to  resolve  gap  widths  on  the  order  of  .2 mm  between  a  dielectric  slab  of 
20c?n  and  a  metal  (perfectly  conducting)  surface  using  an  interrogating  signal  with  a  3mm.  wavelength. 

Future  work  on  this  problem  will  likely  involve  more  efficient  computational  methods  since  currently  the 
inverse  problem  involving  a  20 cm  slab  takes  10  hours.  Further,  more  sophisticated  models  for  describing 
the  polarization  mechanisms  in  non-homogeneous  materials  must  be  developed.  Finally,  in  order  to  take 
scattering  and  non-normally  incident  electromagnetic  signals  into  account,  multi-dimensional  models  will  be 
necessary. 
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